#include "analysis.h"
#include "numutils.h"
#include "globvars.h"
#include "oaffile.h"
#include "gts_tools.h"

#include "loader.h"




CAnalysis::CAnalysis():
  sort_by_id(ps.get<0>()),
  sort_by_snap(ps.get<1>())
#ifdef R_INDEX
,
  sort_by_R(ps.get<2>())
#endif
{
#ifndef DEBUG_ASS

  m_ComVec=new CComVec(string(COM_FILE));
  GenName(fnameApPeaks, TYPE, "peaks_Ap_"); 
  m_ApVec=new CApVec(fnameApPeaks);


 if(!m_ComVec->load_com())
    {	  
      {
	scoped_timer timethis("\t#CAnalysis::SetupCom ");
	TSpecificFunctor<CAnalysis> A(this, &CAnalysis::SetupCom);
	SetNumberOfThreads(8);	
	TFunctor* vAnalysisTable[] = {&A};
	cout<<"START/END: "<<SNAP_START<<" "<< SNAP_END<<endl;	
	parallel_do_RUN(vAnalysisTable[0],SNAP_START, SNAP_END);
	
      }
      
    }

  char buf[2500], snapname[2500];
	
  int isnap=SNAP_START-1;
  bool file_flag=true;
  do{
    if(SNAP_END<=isnap++)break;
    sprintf(snapname, SnapMask, isnap);
    sprintf(buf, "%s/%s", SnapDir, snapname);
    file_flag=FileExists(buf);
    cout<<"#"<<buf<<endl;;
  } 
  while(!file_flag);

  if(!file_flag){cout<<"No Valid imput Files: exiting..."<<endl;stop;};
  SNAP_START=isnap;
  string fname=string(buf);
  cout<<"# reading data from "<<fname<<endl;
  CGetData *data= new CGetData(fname, TYPE);
  ////////////////////////////////////////

  CComPos cpos=GetCom(SNAP_START);

  ////////////////////////////////////////
 
  double RMIN2=RMIN*RMIN, RMAX2=RMAX*RMAX, R2=0.0;
  int  np=data->GetNp();
  cout<<"# Check  particles "<<np<<endl;
  for(int i=0, ic=0; ic<np;ic++)
    {
      //Here we will init   id_sel
      data->pPOS[i]   -= cpos.Pos[0];
      data->pPOS[i+1] -= cpos.Pos[1];
      data->pPOS[i+2] -= cpos.Pos[2];
      i+=3;
    }

 if(OAF_FROM_IDFILE==1){
  CAsciiReader file_ids(string(OAF_ID_FILE), 2);// ID, PosX, PosY
  if(file_ids.IsGood())
    {
       for(TAsciimap::iterator it=file_ids.m_data.begin();it!=file_ids.m_data.end();it++)
	 {
	   cout<<(*it).first<<endl;
	   idu_sel.insert((*it).first);
	 }
    }
    
}else
  {
  for(int i=0, ic=0; ic<np;ic++)
    {
      R2=(data->pPOS[i]*data->pPOS[i]+data->pPOS[i+1]*data->pPOS[i+1]);
      if( R2>=RMIN2 && R2<RMAX2)
	{
	  idu_sel.insert(data->pID[ic]);	  
	}
      i+=3;

    }
  }
  delete data;

  cout<<"# We will keep track for n= "<<GetNTracked()<<
    "  particles "<<endl;
  cout<<"# In range: "<<RMIN<<" -  "<<RMAX<<endl;


#endif
}
  
void CAnalysis::my_print_stats()
{
#ifndef DEBUG_ASS
  if(sort_by_id.size()!=0)
    {
  particle_set::nth_index_iterator<0>::type it_id_b= sort_by_id.begin();
  particle_set::nth_index_iterator<0>::type it_id_e= sort_by_id.end();
  it_id_e--;
	
#ifdef R_INDEX
  particle_set::nth_index_iterator<2>::type it_Rb= sort_by_R.begin();
  particle_set::nth_index_iterator<2>::type it_Re= sort_by_R.end();
  it_Re--;
#endif
	
  snap_unique_set::nth_index_iterator<0>::type it_SUb=snapu.begin();
  snap_unique_set::nth_index_iterator<0>::type it_SUe=snapu.end();
  if(it_SUb!=it_SUe) 
    it_SUe--;
  cout<<"#############################################"<<endl;
  cout<<"Number of Unique Snaps: "<<snapu.size()<<endl;
  cout<<"Number of Unique IDs  : "<<idu.size()<<endl;
  cout<<"Number of loaded IDs  : "<<sort_by_id.size()<<endl;
	
  cout<<"Min  IDs         : "<<(*it_id_b).id<<endl;
  cout<<"Max  IDs         : "<<(*it_id_e).id<<endl;

#ifdef R_INDEX    
  cout<<"Min  R           : "<<(*it_Rb).R<<endl;
  cout<<"Max  R           : "<<(*it_Re).R<<endl;
#endif
	
  cout<<"Min  SNAP        : "<<(*it_SUb)<<endl;
  cout<<"Max  SNAP        : "<<(*it_SUe)<<endl;
    }
  cout<<"#############################################"<<endl;
#endif	
}

int CAnalysis::parallel_hello_data(unsigned short  isnap, unsigned short myThreadID)
{
  {
    boost::mutex::scoped_lock scoped_lock(m_ps_fill_mutex);
    std::cout<<"Hello from:"<<myThreadID<<std::endl;
  }
  return 0;
}
#define LOC_VARS
#define BENCH_FOR_EACH

int CAnalysis::parallel_load_data(unsigned short  isnap, bool append, unsigned short myThreadID)
{  
  if(append)
    {}else;

  int np=0;
  std::string filename;
  GetFileName(filename, isnap);
  CGetData *data= new CGetData(filename, TYPE);
  CComPos cpos=GetCom(isnap);
  np=data->GetNp();
  // data->Translate(cpos.Pos, cpos.Vel);
   for(int i=0, ic=0; ic<np;ic++)
    {
      //Here we will init   id_sel
      data->pPOS[i]   -= cpos.Pos[0];
      data->pPOS[i+1] -= cpos.Pos[1];
      data->pPOS[i+2] -= cpos.Pos[2];
      i+=3;
    }

  ///////////////////////////
  vector<unsigned int> ih;
  particle_set loc_ps;
  id_unique_set loc_idu;
  snap_unique_set loc_snapu;
  time_unique_set loc_timeu;
  ///////////////////////////
  for(int ic=0; ic<np;ic++)
    {
      if( isInSel( data->pID[ic]))
	ih.push_back(ic);
    }
  {
#ifndef LOC_VARS
    boost::mutex::scoped_lock scoped_lock(io_mutex);  
   if(m_verbose) scoped_timer timeme("Inserts: ");
#endif
    ////////////////////////////////////////
   
    char tbuf[5];
    sprintf(tbuf, "%d", isnap);
    string snapio(tbuf);
    int i=0, ic;
    unsigned int isize=ih.size();
    for(unsigned int is=0;is<isize;is++)
      {
	
	  	   
	    
	      ic=ih[is];
	      i=3*ic;
#ifdef LOC_VARS
	      loc_ps.insert(particle(data->pID[ic],isnap,
				 &data->pPOS[i],
				 &data->pVEL[i], 
				 data->pRHO[ic],
				 data->pHSML[ic]));
	    
	      loc_idu.insert(data->pID[ic]);// Unique Sorted Ids by BOOST multi_index		      		      
	      loc_snapu.insert(snapio);// Unique sorted SNAPS by BOOST multi_index
	      loc_timeu.insert(data->GetTime());// Unique sorted TIME by BOOST multi_index
#else	      
	      ps.insert(particle(data->pID[ic],isnap,
				 &data->pPOS[i],
				 &data->pVEL[i], 
				 data->pRHO[ic],
				 data->pHSML[ic]));
	    
	      idu.insert(data->pID[ic]);// Unique Sorted Ids by BOOST multi_index		      		      
	      snapu.insert(snapio);// Unique sorted SNAPS by BOOST multi_index
	      timeu.insert(data->GetTime());// Unique sorted TIME by BOOST multi_index
#endif
	
      }
    {
#ifdef LOC_VARS
      boost::mutex::scoped_lock scoped_lock(io_mutex);
      //     if(m_verbose) 
       scoped_timer timeme("LOCAL Inserts: ");
      
#ifdef BENCH_FOR_EACH
      timeme.SetText(" COPY :");
      std::copy(
		    loc_ps.begin(),
		    loc_ps.end(),
		    std::inserter(ps, ps.end())
		    );
      
      
#else
      ps.insert(loc_ps.begin(), loc_ps.end());
#endif
      idu.insert(loc_idu.begin(), loc_idu.end());
      snapu.insert(loc_snapu.begin(), loc_snapu.end());
      timeu.insert(loc_timeu.begin(), loc_timeu.end());


      //      if(m_verbose) 	
      cout<<"#Done Load: "<<isnap<< " ";


    }
    ////////////////////////////////////////
    delete data;
  }
#endif
  return np;
}
int CAnalysis::load_data(string filename, int TYPE, string snapio)
{  

  int np=0;
  boost::timer t;
  CGetData *data= new CGetData(filename, TYPE);
  {

    std::cout.precision(2);
    std::cout <<"Time using read snapshot: " << t.elapsed();

    int isnap=0;
    //m_type.push_back(TYPE);
    ////////////////////////////////////////
    isnap=atoi(snapio.c_str());
    //  sprintf(buf, "%04d", isnap);
    CComPos cpos=GetCom(isnap);

    np=data->GetNp();
  
    for(int i=0, ic=0; ic<np;ic++)
      {
	if( isInSel( data->pID[ic]))
	  {
	    data->pPOS[i]   -= cpos.Pos[0];
	    data->pPOS[i+1] -= cpos.Pos[1];
	    data->pPOS[i+2] -= cpos.Pos[2];
	    data->pVEL[i]   -= cpos.Vel[0];
	    data->pVEL[i+1] -= cpos.Vel[1];
	    data->pVEL[i+2] -= cpos.Vel[2];
	  
	    ps.insert(particle(data->pID[ic],isnap,
			       &data->pPOS[i],
			       &data->pVEL[i], data->pRHO[ic],data->pHSML[ic]));

	    //ps.insert(particle(data->pID[ic]));


	    idu.insert(data->pID[ic]);// Unique Sorted Ids by BOOST multi_index		      		      
	    snapu.insert(snapio);// Unique sorted SNAPS by BOOST multi_index
	    timeu.insert(data->GetTime());// Unique sorted TIME by BOOST multi_index
	  
	    //	  cout<<ic<<" - "<<np<<endl;
	  }
	i+=3;
      }
    std::cout.precision(2);
    std::cout << "s  insert particle " << t.elapsed() << "s  ";
    ////////////////////////////////////////
    delete data;
  }

  return np;
}

void CAnalysis::OpenOrbitFile(int myID){
  char buf[500];
  sprintf(buf, "%s/OrbitDump_%d", OutputDir,myID);
  ofs.open(buf);
  if(!ofs.good())cout<<"cannot open file: "<<buf<<endl;
  else
    {
      ofs<<"# myID"<<endl;
      ofs<<"# SNAP(1)   ID(2)    X(3)    Y(4)    Z(5)    R(6)"<<endl;
    }
};//prepare Ofs
void CAnalysis::DumpOneOrbit( int inID)
{
#ifdef DEBUG_ASS
  particle_set::nth_index<0>::type& sorted_index=ps.get<0>();
  particle_set::nth_index_iterator<0>::type it1=
    sorted_index.lower_bound(inID);

  OpenOrbitFile(inID);//prepare Ofs
  //boost::archive::text_oarchive oa(ofs);

  while( ofs.good()&& it1 !=sorted_index.end() && ((*it1).id==inID)){
    ////////////////////////////////
    ofs << (*it1);
    ////////////////////////////////
    ++it1;
  }
  ofs.close();
#endif
#ifdef MYDEBUG
  cout<<"#########  END DUMP #######"<<endl;
#endif
}


typedef boost::multi_index::index<CAnalysis::particle_set,id>::type::iterator T_indexID;
struct setSnapLess {
  bool operator( )(const T_indexID p1,
		   const T_indexID p2) {
    return( (*p1).snap < (*p2).snap);
  }
};

void CAnalysis::DumpOrbitsOAF()
{
#ifndef DEBUG_ASS
  /////////////////////////
	if(ps.empty())return;
  string ext=ToString(RMIN)+"_"+ToString(RMAX);
  COAFFile oafFILE(OutputFile+ext);

  //////////////////////////
  snap_unique_set::nth_index_iterator<0>::type it_SUb=snapu.begin();
  snap_unique_set::nth_index_iterator<0>::type it_SUe=snapu.end();
  if(it_SUb!=it_SUe) 
    it_SUe--;
  /////////////////////////////////////////////
	
	
  double T1=0.0;
  double T2=0.0;
	
  unsigned int Ns=snapu.size(), Nsc=0,iskip=0;

  oafFILE.Open('w');
  oafFILE.WriteHeader(idu.size(),string(strOAF_32), T1, T2, 0.005, Ns);

  char *pData;
  //////////////////////////////////
  unsigned int icount=0,nbw=0;
  //////////////////////////////////
  unsigned int timeSize=timeu.size()*sizeof(double);
  pData= new char[timeSize];
  vector<double> Time;
  for(time_unique_set::iterator itt=timeu.begin();itt!=timeu.end();++itt)
    {
      memcpy((void*)&pData[icount], (void*)&(*itt), 8);
      Time.push_back((*itt));
      icount+=8;
    }
	
  oafFILE.WriteTime(pData, timeSize);

  delete [] pData;
  /////////////
  unsigned int DataSize=oafFILE.WriteDataStart(Ns);
  icount=0;
  float var[8];
  pData= new char[4+32*Ns];
  memset(pData, 0, 4+32*Ns);
  set<T_indexID, setSnapLess> pset;
  for(id_unique_set::iterator it=idu.begin();it!=idu.end();++it)
    {
      int myID = (*it); 
      /////////FillOrbitBuffer( myID);
	  T_indexID ic0, ic1;
	  boost::tuples::tie(ic0,ic1)=get<id>(ps).equal_range(myID);

      //particle_set::nth_index<id>::type& sorted_index=ps.get<id>();
      //particle_set::nth_index_iterator<id>::type it1=sorted_index.lower_bound(myID);
	  if(ic0 ==ic1)
	  {
		  cout<<"Strange Error in ID sort"<<endl;
		  stop;
	  }
	
      Nsc=0;
      memmove((void*)&pData[0], (void*)&myID, 4);
      int count=0;
      /* MULTI THREADED THINGS */
      /*We need to sort them by time*/
      pset.clear();		  
      while( ic0 !=ic1 && ((*ic0).id==myID))
	  {
		  pset.insert(ic0);
		  ic0++;
	  }
      
      set<T_indexID, setSnapLess>::iterator itSet=pset.begin();
      while( itSet !=pset.end()){
	////////////////////////////////
	
	var[0]=(*itSet)->Pos[0];
	var[1]=(*itSet)->Pos[1];
	var[2]=(*itSet)->Pos[2];

	var[3]=(*itSet)->Vel[0];
	var[4]=(*itSet)->Vel[1];
	var[5]=(*itSet)->Vel[2];
	
	var[6]=(*itSet)->Rho;	
	var[7]=(*itSet)->Hsml;
	/*std::cout.precision(3);
	  cout<<Time[count++]<<" "<<myID<<" "<<var[0]<<endl;*/
	if(Nsc > Ns*32 )
        {
        cerr<<"ERROR: Nsc > Ns*32 : "<<Nsc<<endl;
        cerr<<"ID: "<<myID<<" Nsc= "<<Nsc<<" we got:"<<Nsc/32<<" should be  Ns= "<<Ns
        <<" SNAP: "<<(*itSet)->snap<<endl;
        cerr<<(*(*itSet));
        //error_flag=1;
        iskip++;
        break;
        }else
	  {
	    
	    memmove((void*)&pData[4+Nsc], 
		    (void*)&var[0], 32);	  
	    Nsc+=32;
	  }
	////////////////////////////////
	++itSet;
      }
      oafFILE.WriteOneParticle(pData, 4+32*Ns);
      nbw+=(4+32*Ns);
      ///////////////////////////////////////
      if(icount%1000 ==0 )//&&VERBOSE_flag)
      {
	cout<<icount<<" ) ";
	cout<<"OAF Dumping the orbit for ID: "<< myID<<endl;
      }
      
      ///////////////////////////////////////
      icount++;
      
    }
  
  oafFILE.WriteDataStop();
  oafFILE.Close();
  delete [] pData;
  std::cout.precision(10);
  cout<<"####################"<<endl;
  cout<<"Dumped n= "<<icount<<" orbits"<<endl;
  cout<<"DATA block total size :"<<nbw/(1024.0*1024.0)<<"Mb"<<endl;
  cout<<"Number of  BAD particles with duplicate IDS: "<<iskip<<endl;
  cout<<"####################"<<endl;

#endif	
    }


  void CAnalysis::DumpOrbits()
    {
#ifdef DEBUG_ASS
      int icount=0;
      for(id_unique_set::iterator it=idu.begin();it!=idu.end();++it)
	{
	  int myID = (*it);
	  string sn=*snapu.begin();
		
	  if(icount%100 ==0)
	    {
	      cout<<icount<<" ) ";
	      cout<<"Dumping the orbit for ID: "<< myID<<endl;
	    }
	  DumpOneOrbit( myID);
	  icount++;
		
	}
      cout<<"####################"<<endl;
      cout<<"Dumped n= "<<icount<<" orbits"<<endl;
#endif
    }

  bool CAnalysis::isInSel(int inID)
    {
#ifndef DEBUG_ASS
      id_unique_set::iterator iidu= idu_sel.find(inID);
      if(iidu==idu_sel.end())
	{
#ifdef VERBOSE
	  std::cout<<"no such ID in the table"<<std::endl;
#endif
	  return false;
	}
#endif
      return true;
    }

  bool CAnalysis::Load(void)
    {
#ifdef DEBUG_ASS
      boost::timer time;
      cout<<"Loading data..."<<endl;
      string dump_file_name=string(OutputFile)+string("_particle_dump");
      std::ifstream ifs(dump_file_name.c_str(),ios::in | ios::binary );
      if(!ifs.good())return false;
      boost::archive::binary_iarchive ia(ifs);
      
      ia>>ps;
      ifs.close();
      for( particle_set::iterator itP=ps.begin();itP!=ps.end();itP++)
	{
	  idu.insert((*itP).id);// Unique Sorted Ids by BOOST multi_index
	  snapu.insert((*itP).snap);// Unique sorted SNAPS by BOOST multi_index
	};
      cout<<"Loaded N= "<<ps.size()<<" particles from dump file."<<endl;
      std::cout << "Time using serialization library: " << time.elapsed() << std::endl;
#endif
      return true;
    }
  void CAnalysis::Save(void )
    {
#ifdef DEBUG_ASS
      boost::timer time;
      cout<<"Strating to dump for warm runs..."<<endl;
      string dump_file_name=string(OutputFile)+string("_particle_dump");
      std::ofstream ofs(dump_file_name.c_str(),ios::out | ios::binary );
      if(!ofs.good()){
	cerr<<"cannot open the dump file"<<endl;return;
      }
      else 
	cout<<"Dump File is ok"<<endl;
      
      boost::archive::binary_oarchive  oa(ofs);
      /* 
     for( particle_set::iterator itP=ps.begin();itP!=ps.end();itP++)
	oa.save_binary(&(*itP), sizeof((*itP)));*/
      oa<<ps;
      /*for( particle_set::iterator itP=ps.begin();itP!=ps.end();itP++)
	{
	  ofs.write((char*)&(*itP).Pos[0], 3*sizeof(float));
	  }*/
      ofs.close();
      cout<<"done dump file  write"<<endl;
      std::cout << "Time using serialization library: " << time.elapsed() << std::endl;
#endif
    }
 //////////////////////////////////////////////////
void  CAnalysis::do_parallel_load_on_SNAPS(int a, int b) 
	  { 
		  for(int isnap=a;isnap<b;isnap++)
		  {
		    this->parallel_load_data(isnap);
		 

		  }
	  };
void  CAnalysis::do_parallel_Profiler_on_SNAP(int a, int b)
	  {
		  for(int isnap=a;isnap<b;isnap+=SNAP_STEP)
		  {
		    this->parallel_profiler_data(isnap);
		  }
	  }
void  CAnalysis::do_parallel_Test_on_SNAP(int a, int b)
	  {
		  for(int isnap=a;isnap<b;isnap+=SNAP_STEP)
		  {
		    this->parallel_Test_data(isnap);
		  }
	  }
 void  CAnalysis::do_parallel_save_on_OAF(int a, int b) 
	  { 
		  for(int id=a;id<b;id++)
		  {
		    this->parallel_load_data(id);		    
		  }
	  };

void CAnalysis::SetupCom(int a, int b)
{
 
  std::string filename;
  float COM[6];
  for(unsigned short isnap=a;isnap<b;isnap++)
    {
      GetFileName(filename, isnap);

      if(FileExists(filename.c_str()))
	{
	  cout<<"  # Before getting data,here file is: "<<filename<<endl;
	  CGetData *data= new CGetData(filename,4);//Warning we are using COM of the Disk Particles
	  cout<<"  # here file is: "<<filename<<endl;
	  data->SetIsnap(isnap);
	  /////////////////////////
	  data->GetComByAFOF(&COM[0], AFOF_MASK);
	  /////////////////////////  
	  delete data;
	  	 	  
	  {
	    boost::mutex::scoped_lock scoped_lock(io_mutex);
	    m_ComVec->PutCom(isnap,CComPos(COM[0], COM[1], COM[2],COM[3], COM[4], COM[5])); 
	    std::cout<<"\t# Loading COM: "<<isnap<<std::endl;	
	  }
	  
	} 
    }


}

CComPos CAnalysis::GetCom(int isnap)
{
   CComPos cpos;   
  if(!m_ComVec->GetCom(isnap,cpos))
    {
      SetupCom(isnap, isnap+1);
      m_ComVec->GetCom(isnap,cpos);
    };

  return cpos;
}

bool CAnalysis::InitData(int isnap, CGetData *&data)
  {
	  std::string filename;
	  float phi=0;
	  GetFileName(filename, isnap);
	  
	  data= new CGetData(filename, TYPE);
	  CComPos cpos=GetCom(isnap);
	  data->Translate(cpos.Pos, cpos.Vel);
	  data->Rotate(phi);

	  return (1 && data->GetNp());
  };





template<class T>
void DumpVector(string filename, std::vector<T> vec, unsigned int n=0,bool verbose=false)
{
  ofstream of(filename.c_str());
  if(!of.good())assert("Error");
  unsigned int i=0, np=min(n,(unsigned int )vec.size());
  for(i=0;i<np;i++)
    {
      if(verbose)cout<<std::setprecision(2)<<std::fixed<<i<<") "<<vec[i]<<endl;
      of<<i<<" "<<vec[i]<<endl;
    }
  of.close();
}
template<class T>
void DumpVector2(string filename,std::vector<T> vecX, std::vector<T> vec, unsigned int n=0,bool verbose=false)
{
  
  ofstream of(filename.c_str());
  if(!of.good())assert("Error");
  unsigned int i=0, np=min(n,(unsigned int )vec.size());
  for(i=0;i<np;i++)
    {
      if(verbose)cout<<std::setprecision(2)<<std::fixed<<i<<") "<<vec[i]<<endl;
      of<<i<<" "<<vecX[i]<<"  "<<vec[i]<<endl;
    }
  of.close();
}
template<class T>
void DumpVector3(string filename,string head, std::vector<T> vecX, 
		 std::vector<T> vec,std::vector<T> vec1,  
		 unsigned int n=0,bool verbose=false)
{
  
  ofstream of(filename.c_str());
  if(!of.good())assert("Error");
  unsigned int i=0, np=min(n,(unsigned int )vec.size());
  of<<head<<endl;
  for(i=0;i<np;i++)
    {
      if(verbose)cout<<std::setprecision(2)<<std::fixed<<i<<") "<<vec[i]<<endl;
      of<<i<<" "<<vecX[i]<<"  "<<vec[i]<<"  "<<vec1[i]<<endl;
    }
  of.close();
}

template<class T>
void DumpVector4(string filename,string head, std::vector<T> vecX, 
		 std::vector<T> vec,  std::vector<T> vec1, std::vector<T> vec2,
		 unsigned int n=0,bool verbose=false)
{
  
  ofstream of(filename.c_str());
  if(!of.good())assert("Error");
  unsigned int i=0, np=min(n,(unsigned int )vec.size());
  of<<head<<endl;
  for(i=0;i<np;i++)
    {
      if(verbose)cout<<std::setprecision(2)<<std::fixed<<i<<") "<<vec[i]<<endl;
      of<<i<<" "<<vecX[i]<<"  "<<vec[i]<<"  "<<vec1[i]<<"  "<<vec2[i]<<endl;
    }
  of.close();
}
template<class T>
void DumpVectorByIndex(string filename, std::vector<T> vec, 
		       vector<unsigned int> idx,
		       unsigned int n=0,bool verbose=false)
{
  ofstream of(filename.c_str());
 
  unsigned int i=0,ic=0, np=min(n,(unsigned int )vec.size());
  for(i=0;i<np;i++)
    {
      ic=idx[i];
      if(verbose)
      cout<<i<<" "<<" "<<ic<<" "<<std::fixed<<endl;
      of<<i<<" "<<" "<<ic<<" "<<endl;
    }
  of.close();
}
enum eMODES
  {
    BY_R,
    BY_RHO,
    BY_ID,
    BY_HSML,
    BY_NONE
  };

template<class T>
bool compR(T v1, T v2)
{
  return (*v1).R < (*v2).R;
}

/*template<class T>
bool compPos(T v1, T v2)
{
  return (*v1).Pos[0] < (*v2).Pos[0];
}*/

 template<class T>
 class compPos
 {
 public:
   compPos(int ind):m_ind(ind){
     if(ind>3 || ind<0)
       {cout<<"Error using template:compPos"<<endl;exit(-13);}
   }
 bool operator()(const T v1, const T v2)const 
   {
     return v1.Pos[m_ind] < v2.Pos[m_ind];
     }
 private:
   int m_ind;
};

template <class T>
class ProxyLess
 {
  T& that;
public:
  ProxyLess(T& f, eMODES mode_=BY_R) : that(f),m_smode(mode_) {}
bool operator()(unsigned int leftID, unsigned int rightID) const 
  {
    
    if(m_smode==BY_R)
      return that.dist(&(that.pPOS[leftID*3])) < that.dist(&(that.pPOS[rightID*3]));

    if(m_smode==BY_RHO)
      return (that.pRHO[leftID] < that.pRHO[rightID]);

    if(m_smode==BY_ID)
      return (that.pID[leftID] < that.pID[rightID]);

    if(m_smode==BY_HSML)
      return that.pHSML[leftID] < that.pHSML[rightID];
    
    return leftID < rightID;
}

 private:
  eMODES m_smode;
};

template <class T>
struct ProxyLog{
  T  operator () ( T tvar){return log(tvar);};
};

//#include "Minimizer.h"
#include "IdTest.h"

void CAnalysis::parallel_profiler_data(int isnap)
{
  

  CGetData *data=NULL;
  if(!InitData(isnap, data))return;
	/////////////////////////////////////////////       
	unsigned int i,ii,np3, np=data->GetNp();
	vector<double> R(np), Rlog(np);	
	np3=np*3;
       

	vector<unsigned int > indexR(np);	

	for(i=0;i<np;i++)
	{
	  R[i]=(data->dist(&(data->pPOS[i*3])));
	  indexR[i]=i;
	}
	
	std::sort( indexR.begin(), 
		   indexR.end(),
		   ProxyLess<CGetData>(*data, BY_R)
		   );
	//std::reverse(indexR.begin(), indexR.end());

	//	DumpVectorByIndex("testi.txt", R,indexR, 10, true);
	string fname;
       
	GenName(fname, isnap, "profile_R_");
	//THistogram<float, float> histo(256);
	//histo.Compute(R);	
	//std::transform(R.begin(), R.end(), Rlog.begin(),ProxyLog<float>());

	//DumpVector2(fname+"Left_",histo.GetLeftContainers(),
	//	    histo.GetNormalizedHistogram(), histo.size(), false);
	//DumpVector2(fname+"Cent_",histo.GetCenterContainers(),
	//	    histo.GetNormalizedHistogram(), histo.size(), false);
	/////////////////////////////////////////////
	//GetAB
	int iR20=count_if(R.begin(), R.end(), 
			  std::bind2nd(std::less_equal<float>(),profRMAX));

	int nshell=static_cast<int>(profRMAX/xStep)+1;

	vector<double> al, al0,al2, bl2, fR;;
	al.resize(nshell);al2.resize(nshell);bl2.resize(nshell);
	al0.resize(nshell);fR.resize(nshell);
	fill(al.begin(),al.end(),0.0 );
	fill(al0.begin(),al0.end(),0.0 );
	fill(al2.begin(),al2.end(),0.0 );
	fill(bl2.begin(),bl2.end(),0.0 );
	fill(fR.begin(),fR.end(),0.0 );
	
	int ishell=0;
	double theta=0.0;
	double massp=1.0;
	//	cout<<iR20<<" "<<R[indexR[0]]<<" "<<R[indexR[iR20]]<<endl;
	for(int i=0;i<iR20;i++)  
	  {
	    ii=indexR[i];
	    ishell=static_cast<int>(R[ii]/xStep);
	    theta=atan2(data->pPOS[ii*3+1], data->pPOS[ii*3]);	  
	    al2[ishell] += massp*cos(fm*theta);
	    bl2[ishell] += massp*sin(fm*theta);
	    al0[ishell] += massp;
	  }

	for(i=0;i<nshell;i++)
	  {
	    al[i]=sqrt(al2[i]*al2[i]+
		       bl2[i]*bl2[i]);
	    if(al0[i]>0.0f)
	      al[i]/=al0[i];
	    fR[i]=xStep*i;
	  }
	double total=accumulate(al.begin(), al.end(), 0.0);
	string  fmFile=fname+string("_Afm")+ToString(fm)+'_';

	DumpVector4(fmFile,"# Rad(kpc)  Al  Alm2  Blm2 ", fR,al,al2, bl2,
		    al.size(), false);


	/////////////////////////////
	//Index sorted by rho
	//	std::sort( indexR.begin(), 
	//	   indexR.end(),
	//	   ProxyLess<CGetData>(*data, BY_RHO)
	//	   );
	// We need sorted by desceding order;
	//	std::reverse(indexR.begin(), indexR.end());
	/////////////////////////////
  //;;;;;;;;;;;;;;;;;;; SET log bins ;;;;;;;;;;;;;;;;;;
	vector<double> hcoef(NCOEF);
	vector<double> hcoef3(NCOEF);
	  
	double dr=log(profRMAX)/(double)(NCOEF);
	for(int i=0;i<NCOEF;i++)
	  {    
     	    hcoef[i]= (exp(dr*(i+1))-1)/(profRMAX);
	    hcoef3[i]= hcoef[i]*hcoef[i]*hcoef[i];
	  }

	const bool NO_MASS_FLAG=true;
	double Rall, R3=profRMAX*profRMAX*profRMAX;
	const double facpi=4.0/3.0*M_PI;
	
	
	vector<double> fRho, fbyDRho;
	fR.resize(NCOEF);
	fRho.resize(NCOEF);fbyDRho.resize(NCOEF);
	fill(fRho.begin(), fRho.end(), 0.0f);
	fill(fbyDRho.begin(), fbyDRho.end(), 0.0f);
	for(int i=0;i<NCOEF;i++)   
	  {
	    if(!NO_MASS_FLAG)
	      {
		Rall=profRMAX*hcoef[i];
		int iR=count_if(R.begin(), R.end(), 
				std::bind2nd(std::less_equal<double>(),Rall));		
		fR[i]+=iR;
	      }
	    
	    /////////////////////
	    Rall=profRMAX*hcoef[i];
	    int i1=count_if(R.begin(), R.end(), 
		     std::bind2nd(std::less_equal<double>(),Rall));
	    Rall=profRMAX*hcoef[i+1];
	    int i2=count_if(R.begin(), R.end(), 
			    std::bind2nd(std::less_equal<double>(),Rall));
	    if(i2 > 0 && i1>0)
	      {
		fRho[i]=(i2-i1)/(facpi*R3*(hcoef3[i+1]-hcoef3[i]));
		/////////////////////
		for(int j=i1;j<i2;j++)
		      fbyDRho[i]+=data->pRHO[indexR[j]];
		fbyDRho[i]/=(i2-i1);
		
	      }
	    /////////////////////	    
	  }
	for(int i=0;i<NCOEF;i++) 
	  {
	    hcoef[i]*=profRMAX;	    
	    
	  }

	if(!NO_MASS_FLAG)
	  DumpVector2(fname+"_cummMlog_",hcoef,fR,
		      fR.size(), false);

	for(int i=0;i<NCOEF-1;i++) 
	  {
	    if(hcoef[i]<R[indexR[0]])
	      {
		//cout<<i<<")Remove:"<<hcoef[i]<<"  "<<R[indexR[0]]<<endl;
		fRho.erase(fRho.begin()+i);
		hcoef.erase(hcoef.begin()+i);
		fbyDRho.erase(fbyDRho.begin()+i);
		i=0;
	      }
	  }
	//	cout<<fRho.size()<<" minr "<<R[indexR[0]]<<endl;
	DumpVector3(fname+"_cummRholog_","# R(kpc)  Rholog  RhologbyDATA ",
		    hcoef,
		    fRho,fbyDRho,
		    fRho.size(), false);
	/////////////////////////////
	// Now we need to test the CSmooth and  CSpline classes
	//	CSpline<double> spline(x,y,NCOEF,"cubic ");
	CSmooth<double> smooth(&fbyDRho, (int)min(fbyDRho.size()*0.2, 2.0));
	vector<double> smVec=smooth.GetSmoothed(6);
	DumpVector3(fname+"_cummRhoSM_","# R(kpc)  RhologbyDATA Smoothed",
		    hcoef,
		    fbyDRho,smVec ,
		    hcoef.size(), false);
	/////////////////////////////////////////////

	

	
	
	safeFreeOne(data);
	
}
/* Write Stellar Rho File into _rho file*/
void  CAnalysis::WriteMeanFile(vector<particle> data,string meanfilename)
	{
	  int type=TYPE;
	  ofstream file;
	  CGadget *tg=new CGadget();
	  file.open(meanfilename.c_str(),ios::binary);
	  unsigned int blsize, i,ninreg=data.size();
	  CGadget::io_header head;
	  float *pHsml=new float[ninreg];
	  float *pRho=new float[ninreg];
	  float *pPos=new float[ninreg*3];
	  float *pVel=new float[ninreg*3];
	  int  *pID=new int[ninreg];
	  file.clear();
	  memset(&head,0,sizeof(head));
	  memset(head.npart,0,sizeof(head.npart));
	  memset(head.npartTotal,0,sizeof(head.npart));
	  head.npart[type]=ninreg;
	  head.npartTotal[type]=ninreg;
	  head.num_files=1;
	  /////////////////////////////////////////////////////
	  blsize=sizeof(head);
	  tg->WriteOneBlock(file,"HEAD", (char*)&head, blsize);
	  //////////////////////////////////////////////////////
	  int k=0;
	  for(i=0;i<ninreg;i++)
	    {
	      pRho[i]=data[i].Rho;
	      pHsml[i]=data[i].Hsml;	
	      pID[i]=data[i].id;
	      pPos[k]=data[i].Pos[0];				
	      pPos[k+1]=data[i].Pos[1];
	      pPos[k+2]=data[i].Pos[2];				
	      pVel[k]=data[i].Vel[0];				
	      pVel[k+1]=data[i].Vel[1];
	      pVel[k+2]=data[i].Vel[2];	

	      k+=3;
	    }
	  //////////////////////////////////////////////////////		
	  blsize=sizeof(float)*ninreg*3;
	  tg->WriteOneBlock(file,string("POS "), (char*)pPos, blsize);
	  tg->WriteOneBlock(file,string("VEL "), (char*)pVel, blsize);	 
	  blsize=sizeof(int)*ninreg;
	  tg->WriteOneBlock(file,string("ID  "), (char*)pID, blsize);
	  blsize=sizeof(float)*ninreg;
	  tg->WriteOneBlock(file,string("RHO "), (char*)pRho, blsize);
	  tg->WriteOneBlock(file,string("HSML"), (char*)pHsml, blsize);		
	  /////////////////////////////////////////////////////
	  file.close();
	  delete [] pPos;		
	  delete [] pVel;
	  delete [] pRho;
	  delete [] pHsml;
	  delete [] pID;
	  delete tg;
	}


//using namespace boost;
using namespace boost::accumulators;

typedef accumulator_set<double, stats<boost::accumulators::tag::tail_quantile<left> > > accumulator_t_left;
typedef accumulator_set<double, stats<boost::accumulators::tag::tail_quantile<right> > > accumulator_t_right;

vector<float> CAnalysis::FillAp(int isnap)
{
  tApPeaks c;
 
  if(m_ApVec->Ap.size()>0)
    return m_ApVec->Ap;
  //////

   BarRot=make_pair(-180.0, 180.0);
   m_OmegaFile="/data1/home/safe/arm2arm/OmegaBarModel7.txt";
   	// Find Apo/pericenter
	
	vector<double> Time, PhiBar;
	typedef map<int, float> MapType;
	double Tcurrent=isnap*dTime, Ts, Te;
	if( ReadOmegaFile(m_OmegaFile))
	  {

	    for( MapType::iterator it = OmegaVec.begin(); it != OmegaVec.end(); ++it ) 
	      {
		Time.push_back((*it).first*dTime);
		PhiBar.push_back( double( (*it).second) );
	
	      }
	    GetOneRot( Time,  BarRot,  Tcurrent, &Ts,  &Te,  PhiBar);
	    int a=int(Ts/dTime);
	    int b=int(Te/dTime);
	    //////// Load stuff in parallel 
	    TSpecificFunctor<CAnalysis> 
	      tLoadSnaps(this, &CAnalysis::do_parallel_load_on_SNAPS);
	    SetNumberOfThreads(8);
	    parallel_do_RUN(&tLoadSnaps, a,b);
	    SetNumberOfThreads(1);
	  
	    /////////////// End load stuff///////////
	    THistogram<double, double> Aphist(NCOEF);
	    Aphist.SetMinSpectrum(0.0);
	    Aphist.SetMaxSpectrum(1.0);
	    vector<float> ApVec;
	    vector<particle> MeanSnapshot, CurrentSnapshot;
	    vector<particle> PhaseTrajectory;
	    for(id_unique_set::iterator it=idu.begin();it!=idu.end();++it)
	      {
		int myID = (*it); 
		T_indexID ic0, ic1;
		boost::tuples::tie(ic0,ic1)=get<id>(ps).equal_range(myID);
		if(ic0 ==ic1)
		  {
		    cout<<"Strange Error in ID sort"<<endl;
		    stop;
		  }
		pair< T_indexID, T_indexID > mma =
		  boost::minmax_element(ic0, ic1);//, compR<T_indexID>);
		double tAp=(*(mma.first)).GetR()/(*(mma.second)).GetR();
		Aphist.Update(tAp);
		//	cout<<tAp<<endl;
		ApVec.push_back(float(tAp));
		///////////////////////////////
		//Get mean values:
		int c=std::distance(ic0,ic1);
		accumulator_t_right acc0( boost::accumulators::tag::tail<right>::cache_size = c );
		accumulator_t_right acc1( boost::accumulators::tag::tail<right>::cache_size = c );

		T_indexID ic=ic0;
		particle meanP(myID);
		while(ic!=ic1)
		  {
		    meanP+=(*ic);
		    float  A[3], r=(*ic).GetR();
		    (*ic).get_velz(&A[0]);
		    acc0(r);
		    acc1(A[0]);
		    if( (*ic).snap == isnap)
		      CurrentSnapshot.push_back((*ic));
		    //if( (*ic).GetR() < 5.0)
		    PhaseTrajectory.push_back((*ic));
		    ic++;		   
		  } 
		meanP *= 1.0/double(std::distance(ic0,ic1));
		MeanSnapshot.push_back(meanP);
		
		double tqP=quantile(acc0, quantile_probability = 0.50 );
		double tqP1=quantile(acc1, quantile_probability = 0.50 );
		//	cout<<myID<<" "<<tqP<<" "<<tqP1<<endl;

		///////////////////////////////
	      }
	    cout<<std::distance(idu.begin(), idu.end())<<endl;

	    WriteMeanFile(MeanSnapshot,"mean_snap_0450");
	    WriteMeanFile(CurrentSnapshot,"curr_snap_0450");
	    WriteMeanFile(PhaseTrajectory,"pahse_trj_snap_0450");



	    m_ApVec->SetAp(&ApVec);
	    m_ApVec->save_ap();
	    vector<double> xAp=Aphist.GetLeftContainers();
	    vector<double> Ap=Aphist.GetHistogramD();
	    //////////////////// CONVOL by savgol ////////
	    vector<double> convAp(Ap.size(), 0.0);
	    convAp=Ap;	    	    	    
	    {//Convolution by given filter
	      int nleft = 5;
	      int nright = nleft;
	      int degree = 0;
	      
	      int savnp=nleft+nright+1;
	      //	      vector<double> filter(savnp+2);
	      //get coeff
	      // Calculate Savitzky Golay kernel;
	      double filter[1000], sum=0.0;
	      savgol(&filter[0],savnp,nleft,nright,0,degree);
	      
	      // for(int co=0,i=1;i<=savnp;i++)
	      //	 {
	      //	 std::cout.precision(4);
	      //	 cout<<co++<<" "<<filter[i]<<endl;
	      //	 sum+=filter[i];
	      //	 }
	      //cout<<"sum="<<sum<<endl;
	      
	      int Nlen=convAp.size();

	      int M=savnp;
	      int N=Nlen;
	      int i, j;
	      double cmp;
	      vector<double> data(Nlen+1);
	      copy(Ap.begin(), Ap.end(),&data[1]);

	      for (i=1;i<=Nlen;i++) 
		{
		  cmp=0.0;
		  for (j=1;j<=M/2;j++) 
		    {
		      //     if((i-j)>0 && (i-j)<Nlen)
			{
			  cmp += data[((i-j-1+N) % N)+1]*filter[j+1];
			  cmp += data[((i+j-1) % N)+1]*filter[M-j+1];
			}
			//else{cmp=0;}
		    }
		  cmp += data[i]*filter[1];
		  convAp[i-1]=cmp;	
		}
	    }
	   
	    ////////////////////////////////////////////////
	    CSmooth<double> smoothAp(&convAp, 20);
	    vector<double> smAp=smoothAp.GetSmoothed(15);
	    string tfname;
	    GenName(tfname, isnap, "_ApHist_");
	    cout<<"#### dumpdata: "<<tfname<<endl;
	    DumpVector3(tfname,"# ind  xAp  Ap ",
			xAp,
			smAp,convAp,
			xAp.size(), false);
	    
	    
	    /////////////////////////////////////////////
	    cout<<"# Report for SNAP = "<<isnap<<endl;
	    CPeaks<double> peaks(&xAp, &smAp);
	    peaks.BrutPeaks();
	    peaks.show();
	    /////////////////////////////////////////////
	    /////////////////////////////////////////////
	    vector<pair<double, double> >::iterator itP=peaks.final_peaks.begin();
	    while(itP!=peaks.final_peaks.end())
	      {
		c.push_back(tPeaks((*itP).first,(*itP).second, 1)) ;
		itP++;
	      }
	    c[1].peak=0;
	    {
	      boost::mutex::scoped_lock scoped_lock(io_mutex);
	      m_ApVec->AddToApPeaks(isnap, c);	   
	    }



	    return ApVec;
	  }

}
#ifdef HAVE_GTS
#include "gts_tools.h"
#endif
void CAnalysis::parallel_Test_data(int isnap)
{ 

#ifdef HAVE_GTS
  CGTSRead gts_Object("/net/direct/dnas01/arm2arm/DATA/LIA/SPHBAR/NFW/MODEL7/NEW/SURF/snap_tr_0450.gts");
  gts_Object.Stats();
  // gts_Object.ShowVertexes();
#endif

  CGetData *data=NULL;
  if(!InitData(isnap, data))
    {cout<<"#cannot find file:"<<isnap<<endl;
    return;
    }
  /////////////////////////////////////////////       
  unsigned int i,ii,np3, np=data->GetNp();
  float * pPOS=&data->pPOS[0];
  int * pID=&data->pID[0];
#ifdef HAVE_GTS
  std::vector<int> in_surface = gts_Object.GetPointsInside(&pPOS[0], np);
  cout<<"# number of points inside the surface="<<in_surface.size()<<" out of Np="<<np<<"\n\n\n"<<endl;
#endif 
 for(int i=0;i<np;i++)
    id_set.insert(
		  id_locator(i,pID[i], false, false)
		  );
#ifdef HAVE_GTS
  int is;
  BOOST_FOREACH( is, in_surface)
    {      
      id_set.modify(id_locator(is,pID[is], true, false), SURF_CASE);
    }
#endif
  ////
  vector<float> ApVec = FillAp(isnap);  
  tApPeaks PeakI;
  if(! m_ApVec->GetAp(isnap, PeakI)) {
    cout<<"Sonething odd here"<<endl;
    exit(0);
  }
  BOOST_FOREACH(tPeaks p , PeakI)
    cout<<p;
  if(PeakI.size()!=3){
    cout<<"Sonething odd here, more than 3 peaks in Ap distribution.\nProbably BUG in CPeaks class" <<endl;
    //  exit(0);
  };  
  cout<<"# Coolest peak=";
  cout<<PeakI[1].pos[0]<<endl;

 
 for(int i=0;i<np;i++)
   {   
 
     if(ApVec[i]<PeakI[1].pos[0])
       {

       id_set.modify(
		     id_locator(i,pID[i], false, true), AP_CASE
		     );
       }
       
   }
 string tfname;
 GenName(tfname, isnap, "_ApID_");
 cout<<"# try to dump data to file:\n #  "<<tfname<<endl;
 ofstream fAp(tfname.c_str());
 i=0;
 for(id_unique_set::iterator it=idu.begin();it!=idu.end();++it)
   {
     
     fAp<<ApVec[i]<<" "<<(*it)<<endl;
     i++;
   }
 delete data;

}
